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Abstract 

We calculate the leading mass-suppressed two-loop contribution, proportional to 
Og/(m g r 2 ), to the potential of a heavy quark-antiquark system. This contribution 
originates in the nonabelian nature of quantum chromodynamics (QCD) and has 
no analogue in quantum electrodynamics (QED). For this purpose, we elaborate 
a technique based on the implementation of the threshold expansion within the 
effective theory of nonrelativistic QCD (NRQCD). We discuss phenomenological 
implications of our result for heavy-quarkonium physics. We also confirm a previous 
result for the two- loop 0{a 2 s ) correction to the static potential. 
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Since the pioneering work by Susskind way back in 1977 the static potential 
between a heavy quark and its antiquark has attracted the attention of many authors 
0J3|^|,||,[],[7||| . The structure of the heavy- quark- antiquark potential is crucial both for 
understanding fundamental properties of QCD, such as confinement, and for describing 
the rich phenomenology of heavy quarkonia. While the potential, of course, cannot be 
found in a closed form, it can be systematically computed as an expansion in the strong- 
coupling constant a s and the inverse of the heavy-quark mass m q or, equivalently, in the 
heavy-quark velocity v. The leading-order (LO) approximation is given by the simple 
generalization of the Coulomb potential of QED. The purely relativistic 0(v 2 ) correction 
is given, up to the color factor, by the standard Breit Hamiltonian. The one-loop 0(a s ) 
correction has been known for a long time 0, while the two-loop 0(a 2 ) correction has 
been computed only recently J?]|§. The leading and subleading mass- suppressed one- 
loop 0{ya s ) and 0{v 2 a s ) corrections have been found in Refs. f|||||. In this Letter, we 
compute the two-loop leading mass-suppressed 0{va 2 ) term. 

Let us take a closer look at this specific class of corrections. Owing to the gluon self- 
coupling, the QCD potential differs from its abelian counterpart in two principal respects, 
namely by the running of the coupling constant Jl| and the presence of generic corrections 
non-analytic in v 2 . In general, terms non-analytic in v 2 can also appear in QED as a result 
of iterations of operators that vanish on shell, in connection with the Coulomb singularity, 
which brings in a factor oil/ v. One can avoid such terms by including off-shell operators 
in the Hamiltonian, which is the standard practice in QED bound-state calculations. The 
non-analytic terms that we consider here are a consequence of the nonabelian structure of 
QCD, and they are not necessarily related to the Coulomb singularity. They first appear 
at one loop as the well-known nonabelian a 2 /(m q r 2 ) potential |4j]. In this Letter, we 
take the next step and calculate this potential at two loops. Our main result is given 
by Eqs. ©-(0). In order to derive it, we elaborate a technique based on the effective 
theory of NRQCD flO| implemented with the threshold expansion offering the possibility 
to systematically expand Feynman diagrams near threshold ||. We expect this approach 
to be very useful for solving a wide class of problems of nonrelativistic dynamics. 

The perturbative QCD potential is of special interest, since it is the basic ingredient for 
the theory of heavy quarkonia. In fact, the study of nonrelativistic heavy quark- antiquark 
systems |T]J and its applications to bottomonium [I2j and top-antitop [[L3|] physics rely 
entirely on first principles of QCD. In principle, these systems allow for a perturbative 
treatment, the nonperturbative effects fl4| being well under control and model dependence 
being immaterial. This makes heavy quark-antiquark systems to be an ideal laboratory 
to determine the fundamental parameters of QCD, such as a s and m q . Recently, essential 
progress has been made in the theoretical investigation of the nonrelativistic heavy-quark 
threshold dynamics based on the effective-theory approach |IIJ. The analytical results 
for the main parameters of the nonrelativistic heavy quark-antiquark system are now 
available through next-to-next-to- leading order (NNLO) in a s and v [|T5| , p^| , ^7l , p^| , p^| , and 
they have been applied to bottomonium PE|JTP| , |2T| and top-antitop [f2TJ phenomenology. 
Some specific classes of next-to-next-to-next-to-leading order (N 3 LO) corrections have 
also been investigated p2.E3.l33j]; for a brief review, see Ref. |25||. These corrections have 
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turned out to be so sizeable that it appears to be indispensable to gain full control over 
this order, both in regard of phenomenological applications and in order to understand 
the structure and the peculiarities of the nonrelativistic expansion. This Letter provides 
a major step in this direction. 

In the following, we briefly outline the main features of our method, which is based 
on the combination of the concepts of threshold expansion || and effective theory |K| . 
The effective theory is designed to separate the scales and to perform the expansion in v 
at the level of the Lagrangian. Let us recall that the dynamics of a nonrelativistic quark- 
antiquark pair involves four different scales, or regions ||: (i) the hard scale (energy and 
three-momentum scale like m q ); (ii) the soft scale (energy and three-momentum scale like 
(iii) the potential scale (energy scales like m q v 2 , while three- momentum scales like 
and (iv) the ultrasoft scale (energy and three-momentum scale like m q v 2 ). The 



m q v 



m q v 



ultrasoft scale is only relevant for gluons. NRQCD |T0| is obtained by integrating out 
the hard modes. Subsequently integrating out the soft modes and the potential gluons 
results in the effective theory of potential NRQCD (pNRQCD) [p6| , which contains po- 
tential quarks and ultrasoft gluons as active particles. The dynamics of a nonrelativistic 
quark- ant iquark pair in pNRQCD is governed by the effective Schrodinger equation and 
by its multipole interactions with the ultrasoft gluons. The corrections from the modes 
integrated out are contained in the higher-dimensional operators of the nonrelativistic 
Hamiltonian, corresponding to an expansion in v, and in the Wilson coefficients, which 
are expanded in a s . Both NRQCD and pNRQCD have specific Feynman rules, which 
can be used for the computation of a systematic perturbative expansion. However, this 
problem is complicated because the expansion of the Lagrangian corresponds to a par- 
ticular subspace of the total phase space. Thus, in a perturbative calculation within the 
effective theory, one has to formally impose some restrictions on the allowed values of the 
virtual momenta. Explicitly separating the phase space introduces additional scales to 
the problem, such as momentum cutoffs, and makes the approach much less transparent. 
A much more efficient and elegant method is based on the expansion by regions PJ27|, 
which is a systematical method to expand Feynman diagrams in any limit of momenta 
and masses. It consists of the following steps: (i) consider various regions of a loop mo- 
mentum k and expand, in every region, the integrand in Taylor series with respect to the 
parameters that are considered small there; (ii) integrate the expanded integrand over the 
whole integration domain of the loop momenta; (iii) put to zero any scaleless integral. In 
step (ii), dimensional regularization, with d = 4 — 2e space-time dimensions, is used to 
handle the divergences. In the case of the threshold expansion in v, one has to deal with 
the four regions and their scaling rules enumerated above. 

In principle, the threshold expansion has to be applied to the Feynman diagrams of 
full QCD. However, as we are only interested in the soft and potential contributions, 
it is possible to apply step (i) to the diagrams constructed from the NRQCD Feynman 
rules. Equivalently, the Lagrangian of the effective theory can be employed for a per- 
turbative calculation without explicit restrictions on the virtual momenta if dimensional 
regularization is used and the formal expressions derived from the Feynman rules of the 
effective theory are understood in the sense of the threshold expansion. The last remark 
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is crucial because, in general, the naive use of the effective-theory Feynman rules and 
dimensional regularization leads to an incorrect result. Although this effect has not yet 



shown up in the QED bound-state calculations performed along these lines so far p8|j21 
it is important for our analysis. 

As for the heavy-quark-antiquark potential, the contribution from the hard regions 
is analytic in m q and starts, for dimensional reasons, at order 1/m 2 . By definition, 
the effect of the ultrasoft modes should not be included in the potential, so that we 
are left with the contributions of the soft and potential regions. In the effective-theory 
language, we study the reduction from NRQCD to pNRQCD and compute the effect of 
the soft and potential modes being integrated out. For this purpose, it is enough to use 
the Lagrangian of NRQCD. Apart of the standard LO terms involving the gluon, ghost, 
quark, and antiquark fields in the Lagrangian, we need to also include the l/m g -suppressed 
terms originating from the covariant-derivative operator D 2 / (2m g ) acting on the quark 
and antiquark fields. Note that this operator also includes the quark kinetic-energy term 
k 2 /(2m q ). It can be can be either treated as a perturbation if k is soft or kept in the 
nonrelativistic quark propagator 

S(k) ' 



k — k j (2m q ) + is 

if k is potential. This issue is discussed below in more detail. 

By using the NRQCD Feynman rules, we immediately recover the well-known one-loop 
result, 

where Cp = (N 2 — 1)/(2N) and Ca = N are the eigenvalues of the quadratic Casimir 
operators of the fundamental and adjoint representations of the SU(N) color group, re- 
spectively. In coordinate space, Eq. (||) takes the form 

The abelian part enters the one-loop expression because we use the on-shell potential. In 
fact, we shall need it for the renormalization of our two-loop result. In the standard QED 
analysis using the Coulomb gauge, it is removed by introducing in turn the tree-level 
operator 

v- = =*pl^H, (4) 

where p and p' are the relative three-momenta of the quark and antiquark, respectively, 
and q = p' — p. This is a purely off-shell operator, which vanishes for on-shell quarks, 
when p 2 = p' 2 = 2m q E, where E is the quark energy counted from the threshold. By 
using the Coulomb equation of motion, it is straightforward to check that the matrix 
elements of the abelian part of Eq. (|||) and of Eq. (Q) between Coulomb states are the 
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same. The use of off-shell operators is advantageous in QED because it allows one to 
reduce the number of loops by means of the Coulomb equation of motion, as may be seen 
by comparing Eqs. (Q) and (|J). However, we use the on-shell formulation and the general 
covariant gauge, which is more suitable for multiloop QCD calculations. 

The calculation of the nonabelian part of Eq. (Q) can be greatly simplified by observing 
that the propagator of Eq. ([]]) can be expanded in 1/ m q not only in the soft region, but 
also in the potential one, which is related to the contribution of its pole to the integral over 
k . This is possible because this part is free of Coulomb singularities. In fact, performing 
the expansion, one recovers the familiar generalized functions of ko, S^ n \ko). By contrast, 
the abelian part suffers from the Coulomb singularity corresponding to the iterations of 
the off-shell operator of Eq. (|4]) and the Coulomb potential. Note that the iterations of 
the operators included in the effective Hamiltonian are taken into account by solving the 
Schrodinger equation perturbatively around the Coulomb solution. However, since we do 
not include the off-shell operator of Eq. in the effective Hamiltonian, we must include 
the result of its iterations directly in the potential. After expanding the quark propagator, 
this produces ill-defined products like l/[{k + ie) m (k — ie) n ], which must be avoided by 
resorting to the unexpanded version of Eq. (P in the potential region. 

The structure of the expansion remains intact at two loops, and our final result reads 



V!%,(9) 



na^CF 
m q \q\ 



(C F - 2C A )(3 



4(2CV + C A )C A 



In 



+ b 2 



(5) 



where 



. (HI + Un2) C 2 A + -C A T F n h 
V 36 3 J A 36 



(6) 



j3 = HC A /?> — 4T F ni/2> is the one-loop coefficient of the QCD (3 function, Tp = 1/2 is the 
index of the fundamental representation, and ni is the number of light-quark flavors. For 
iV = 3 and n x = 3, 4, and 5, we have b 2 « -20.836, -18.943, and -17.049, respectively. 
In Eq. (Bp, we omitted the purely Abelian part, since the QED potential to this order is 
known in the off-shell form. The result evaluated in dimensional regularization may be 

m 



found, for example, in Ref. 



In coordinate space, Eq. (El) becomes 



(2) 
/m q 



a^C F 
2irm q r 2 



(C F - 2C A )P 



4(2CV + C A )C A 



In (/iV) + b 2 



(7) 



where \x = e lE fi, with being Euler's constant. In our two-loop analysis, we again used 
the expanded form of Eq. ([!]) in the calculation of the maximal nonabelian structures, 
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CfC\ and CfC 'aTfUi, while it was necessary to keep Eq. ([!]) unexpanded for the treatment 
of the Coulomb singularities associated with the C f Ca and CpTprii terms. 

We performed a number of nontrivial checks for our analysis, (i) We worked in the 
general covariant gauge and verified that the gauge parameter cancels in our final result, 
(ii) The two-loop expression from which Eq. (||) is obtained contains both ultraviolet 
(UV) and infrared (IR) divergences. The UV ones were removed in Eq. (^) by the renor- 
malization of a s in the one- loop result of Eq. (0), which we performed in the modified 
minimal-subtraction (MS) scheme. By the same token, the renormalization group loga- 
rithms proportional to (3 in Eq. compensate the /i dependence of Eq. (§). On the 
other hand, the IR divergences are canceled by the UV ones of the ultrasoft contribution 



22| leaving a finite, //-independent result for the spectrum, (iii) To test our program we 



also recalculated the two-loop correction to the static heavy-quark-antiquark potential 
and found agreement with Ref. ||. This is useful in its own right, since Ref. |§ disagrees 
with the original result of Ref. J7| . 

In Eq. (^|), the IR divergences were subtracted according to the MS prescription. 
For consistency, the same prescription must be used for the calculation of the ultra- 



soft contribution. We should mention that, in Ref. f22fl , a slightly different regulariza- 
tion/subtraction scheme was used. There, only the divergent ultrasoft integrals were 
dimensionally regularized, while the remaining convergent integrals over the potential re- 
gion were solved in four dimensions, despite overall factors containing poles in e. This 
means that an additional matching term must be added to the result of Ref. |22| so as to 
convert it into the conventional MS scheme used in the present paper. In particular, this 
term contains logarithms that exactly cancel one half of the IR logarithms in Eq. (^). The 
remaining IR dependence is canceled by the ultrasoft logarithms, involving the scale E 
instead of |q|, that result in the logarithmically enhanced corrections to the bound-state 
parameters [ [23] . 

Our result is very relevant for heavy-quarkonium spectroscopy. In fact, the complete 
N 3 LO analysis of the spectrum includes three basic ingredients: (i) the leading retardation 
effect; (ii) the multiple iterations of the low-order corrections to the potential; and (iii) 
the N 3 LO potential, which includes (D(v 2 a s ), 0(va 2 ), and 0(a^) terms. The retardation 
effect related to the radiation and absorption of ultrasoft gluons, which results in the QCD 
analogue of the well-known Bethe logarithms, was analyzed in Ref. |22| . The corrections to 
the spectrum due to the multiple iterations of the low-order corrections to the potential can 
be obtained using the technique of Ref. fl8fl . So far, only the one-loop 0{v 2 a s ) part of the 
N 3 LO potential has been available ||,[5]||. Now that also the 0{va 2 s ) part is known, only 
the three-loop C(a 3 ) correction to the static potential remains to be computed. Although 
the latter analysis is tedious from the technical point of view, it is rather straightforward 
from the conceptual one. We believe that the calculation of the two-loop 0(va 2 ) term 
performed here is probably most involved conceptually in the sense that the full power of 
modern technology in effective field theory is required. 

Our result is also very important in view of the recent analysis of Ref. pi| . There, 
it was conjectured that the perturbation theory for the normalization of the ti produc- 
tion cross section at threshold can be essentially improved by reordering the corrections 
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and converting the standard perturbative expansion in powers of a s and v around the 
Coulomb solution to a logarithmic one, which sums up the infinite number of logarithmi- 
cally enhanced terms in each order by solving the nonrelativistic renormalization group 
equation. This conjecture is supported by an analysis involving the known part of the 
NNLO anomalous dimensions, which is, however, incomplete. The remaining part of 
the next-to-next-to-leading logarithmic (NNLL) corrections can be obtained using an ap- 
proach |29| developed for the calculation of the subleading logarithmic corrections in the 
spectrum of positronium bound states. The most nontrivial ingredient necessary to com- 
plete the calculation of the NNLO anomalous dimensions is the two-loop result given in 
the present paper. Note that the QED analysis [29| implies that the missed terms are 
likely to be important. In fact, their abelian counterparts are responsible for the third- 
order subleading logarithmic corrections to the positronium spectrum ||29|| . If the complete 
analysis supported the conjecture of Ref. p4|, this would be crucial for the precision study 



of Higgs-boson-induced effects in tt threshold production. 
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